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Abstract. We introduce a model inspired from statistical physics that is shown to display 
flexible short-range spatial correlations which are potentially useful in geostatistical modeling. 
In particular, we consider a suitably modified planar rotator or XY model, traditionally used 
for modeling continuous spin systems in magnetism, and we demonstrate that it can capture 
spatial correlations typically present in geostatistical data. The empirical study of the spin 
configurations produced by Monte Carlo simulations at various temperatures and stages in the 
nonequilibrium regime shows that their spatial variability can be modeled by the flexible class 
of Matern covariance functions. The correlation range and the smoothness of these functions 
vary significantly in the parameter space that consists of the temperature and the simulation 
time. We briefly discuss the potential of the model for efficient and automatic prediction of 
spatial data with short-range correlations, such as commonly encountered in geophysical and 
environmental applications. 


1. Introduction 

Modern remote sensing techniques allow recording enormous amounts of spatial data, including 
geographical, natural resources, land use, and environmental remote sensing images, which need 
efficient (preferably real-time) processing [IJI2J. Such processing involves the reconstruction of 
missing data that often occur due to different reasons, such as equipment malfunctions and 
gaps in the coverage of the targeted area that appear as a result of restricted satellite paths or 
bad weather conditions |31 HI El El E] . However, such massive data sets can not be efficiently 
handled by standard geostatistical methods, such as kriging [5]. The main drawbacks of such 
methods are high computational complexity, difficulties to automatize the algorithms to work 
without subjective user inputs (selection of variogram and search radius for interpolation) and 
often the necessity of some data pre-processing (e.g., lognormal transformation) if the data does 
not comply with the Gaussianity requirement [3]. Thus, there is a need to develop new spatial 
prediction techniques that overcome these shortcomings mm- Recently, we have proposed 
efficient spatial classification methods based on models inspired from statistical physics, in 
particular discrete spin Ising, Potts, and clock models, employing a heuristic “energy matching” 

principle mm- 

In the present study we extend the idea of using spin models for spatial prediction purposes 
and introduce a new method that is based on a continuous spin planar rotator model. Even 


though the above mentioned nonparametric discrete-spin-based classification models have been 
shown to be rather competitive, their performance was based on the assumption of the existence 
of relevant correlations at some unknown parameters and the prediction results were discrete 
values even for continuous data. As we empirically demonstrate, the modified planar rotator 
model inherently displays flexible short-range spatial correlations that vary significantly over the 
model’s parameter space and could be used for spatial interpolation of continuous geostatistical 
data. 

2. Model and simulations 

2.1. Standard planar rotator model 

The Hamiltonian of the standard two-dimensional planar rotator model with nearest-neighbor 
interactions on a square lattice is defined as 

h = -j^sisj = -jyy° s ((/>j - cf)j), (i) 
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where Si — (cos <^, sin <^) is a continuous spin on i-th lattice site, represented by a two- 
dimensional unit vector, E [0, 2n] is an angle associated with the spin J is an exchange 
interaction parameter and (i,j) denotes the sum over nearest neighbors. The Mermin-Wagner 
theorem m prevents any long-range ordering at finite temperatures in such a model, which has, 
however, been intensively studied in connection with quasi long-range ordering (the so called 
Kosterlitz-Thouless phase) that appears at low temperatures |15. [16]. The phase transition 
is produced by the unbinding of vortex-antivortex pairs at the Kosterlitz-Thouless critical 
temperature Tkt> below which all spins are almost aligned even though true long-range order is 
destroyed by spin fluctuations. The quasi long-range ordering for T < Tkt is characterized by 
the power-law decaying correlation function 
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C(h;n) = {s(h)s( 0)) = ( z ) , (2) 

where h is the spin separation distance, L is the linear lattice size, and 77 (T) = T/2nJ is the 
temperature-dependent exponent. 

2.2. Modified planar rotator model 

Even though the algebraic correlations with the tunable exponent 77 (T) could be relevant for 
modeling of some spatial processes, there are some deficiencies that prevent the straightforward 
use of the standard planar rotator model for such purposes. First of all, one needs to define a 
proper mapping between the spin values and the geostatistical values. Geostatistical data are 
often positively correlated reflecting a degree of spatial continuity, which means that neighbors 
with similar values are more likely (i.e., have lower energy) than those with different values. 
Apparently, this condition is not fulfilled in the standard planar rotator model, described by 
the Hamiltonian (pQ); the latter, for example, assigns equal energies to a pair of neighbors whose 
angles differ by 9 as to a pair whose angles differ by 27 t + 9. This degeneracy could be fixed by 
defining the energy functional so that it monotonically increases with the turn angle between 
the neighboring spins on the entire interval [0, 2 tt] . The latter can be achieved by modifying the 
Hamiltonian (pQ ) to take the following form of the modified planar rotator (MPR) model 

h' = - J Yi cos te(^ _ ( 3 ) 
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where 0<g<l/2isa modification factor. In the following we will use a fixed value of q = 1/2. 
The difference between the nature of spatial correlations in the standard and modified models 



Figure 1 . Snapshots of spin configurations in the low-temperature regime (T = 10 4 ) for (a) 
standard and (b) modified planar rotator model, with periodic boundary conditions. 


is apparent from snapshots of spin configurations (turn angles) in the low-temperature region, 
shown in Fig. [TJ While the presence of the vortex-antivortex pairs shows up in Fig. |l(a)| in 
the form of sharp boundaries between the domains of similarly oriented spins, the variation of 
the spin values in the MPR model, shown in Fig. |l(b)| is rather smooth, as it is typical for 
geostatistical data. 

2.3. Monte Carlo simulations 

We use Monte Carlo (MC) simulations with the Metropolis update rule and vectorized 
checkerboard algorithm. The results are presented for a square lattice with the size L x L. 
We chose L — 128 as a compromise value between the computational resources needed for MC 
simulations and calculations of the empirical variograms (see below) on one side and the effort 
to secure ergodic conditions on the other side. In typical simulations of magnetic systems one 
would like to suppress boundary effects by employing periodic boundary conditions (see the 
snapshots in Fig. [Tj) . However, these are not appropriate for simulation of geostatistical data, 
since normally there is no reason to assume the data on the opposite boundaries are correlated. 
Free boundary conditions are not appropriate either since the spatial process is not necessarily 
confined within the domain boundaries. Instead, we consider it reasonable to apply boundary 
conditions that assume the smooth continuation of the spatial process beyond the borders. This 
assumption is implemented by inserting additional nodes that decorate the lattice and requiring 
that these additional points have the same values as their nearest neighbors inside the lattice. 
Therefore, if Sij is a spin in the i -th row and j -th column of the lattice, i, j = 1,... L, then 
s i,L+i — s L+i,j = slj , = &i,i and sqj = sij, where the rows and columns with indices 

0 and L + 1 respectively are formed by auxiliary nodes around the lattice added to deal with 
the boundary effects. 

2.4- Modeling spatial variability 

The MC simulations produce spatially correlated spin realizations which can be represented 
by their turn angles fa. It turns out that there is great flexibility in the spatial correlations as we 
move in the temperature and simulation time parameter space. Therefore, in order to model the 
spatial variability (the local variogram) of the simulated spin values we use a very flexible Matern 
covariance model. In addition to the trivial parameters of the standard models (Gaussian, 
exponential, spherical, etc.), controlling the variance, cr 2 , and the characteristic covariance 






















length, £, the Matern model involves one more parameter, z/, controlling the smoothness of the 
spatial process. Therefore, it can be considered appropriate for those geostatistical applications 
in which controlling the smoothness is important. Furthermore, the Matern model includes the 
Gaussian and the exponential models as special cases for v —>► oo and v — 0.5, respectively, as 
well as several other bounded models £3- Due to its flexibility, it has been used in various 
areas including pedology PH HUE], hydrology 123 ED, topography ED, health modeling [ 23 ], 
meteorology [23] and environmental modeling [[25]. The Matern correlation function has the 
general form 
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where 6' — (£, v) and K v is the modified Bessel function of order v. Then the corresponding 
variogram function 7 (/i; 0 ), which is typically used in geostatistics, is related to the correlation 
function as 

l(h; 9) = al + a 2 [l - C(h; 6’)] , (5) 

where a 2 is the random field variance, is the nugget variance that corresponds to uncorrelated 
fluctuations and 6 = (cr n ,<7, £, v) is a vector of the complete model parameters. 

Having generated realizations of spin values from MC simulations, we can assess their spatial 
variability by calculating the experimental (i.e., sample-based) variogram as 


N(h) 

7(0 = 2N(h) S “ ^ Xi + O] 2 , (6) 

where N(h) is the number of pairs of spins separated by the distance h. 

There are several methods of fitting the model to the experimental variogram. In the present 
study we used the weighted least squares estimator (hereafter WLS) that was reported to give 
the best overall results m ■ The WLS estimator is based on minimizing the weighted sum of 
squared residuals (objective function), which is defined as follows: 


0 W=;P) 

where 7 (h^) is the experimental and 7 (h^\6) the model variogram. The summation runs over 
the k = 1 ,..., fc max lags containing N(hk) pairs of points, where /i(fc max ) is less than one half of 
the maximum lag in the sample. 


3. Results and discussion 

In order to empirically study spatial correlations in the realizations of spin values generated by 
the modified planar rotator model we run MC simulations at various temperatures T (in units 
J/ks , where ks is Boltzmann constant). Each simulation starts from an initial configuration of 
spatially uncorrelated random spin angles uniformly distributed in the interval [ 0 , 2i r] and the 
snapshots are collected after r MC sweeps. Then, using Eq. ([ 6 1) the experimental variogram 
is calculated and subsequently fitted to the appropriate model by the WLS method defined by 
Eq. ([7j). The results obtained at relatively higher temperatures are demonstrated in Fig. [2j 
for T = 10 -2 and r = 10, 100, 1000 and 10000. Soon after the start of the simulation one can 
observe the development of spatial correlations both in the snapshots in the form of nucleation of 
small spin domains with similar values as well as in the small-scale behavior of the experimental 
variograms. The fitted parameters £ = 1.62 and v — 0.62 for r = 10 in Fig. 2(a) [ indicate that the 
realization is quite rough with a rather small characteristic length. As the relaxation proceeds 
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Figure 2. Experimental and fitted variograms of spin realizations —shown in the inset— 
obtained at temperature T = 10 -2 and at various simulation times r equal to (a) 10, (b) 100, 
(c) 10 3 and (d) 10 4 MC sweeps. Mate(0) denotes the Matern model for the inferred parameter 
values 6 = (<r n , <r, £, v). 


the variance decreases and both the characteristic length and the smoothness parameter increase, 
as shown in Figs. |2(bj| and [2(e)] , for r = 100 and 1000, respectively. Nevertheless, as equilibrium 
is approached the realizations become rougher again, whereas the characteristic length further 
increases and saturates only in equilibrium to a temperature-dependent value. The estimated 
characteristic length £ = 56.51 for the equilibrium configuration at r — 10000 in Fig. |2(d)~] implies 
that the samples on the 128 x 128 lattice suffer from lack of ergodicity, which is also reflected 
in the trend displayed by the corresponding experimental variogram. 

As the temperature is decreased the trend of building up correlations in the equilibration 
process is further enhanced. Since thermal fluctuations are gradually suppressed the realizations 
become smoother. This tendency is apparent by visual comparison of the snapshots obtained at 
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Figure 3. The same as in Fig. [2] for T = 10 4 . Expo($) and Gaus($) denote respectively the 
exponential and Gaussian models for the inferred parameter values $ = (<r n ,<r,£). Note that in 
Fig. |3(a)| the fitted Matern and exponential models coincide as do in Fig. |3(d)| the Matern and 
the Gaussian models. 


T — 10 -4 , shown in Fig.[3]with those generated in the same stages of the relaxation at T — 10 -2 , 
shown in Fig. [2j The difference between the inferred values of the smoothness parameter for the 
respective cases becomes evident particularly at later stages (r = 10 3 and 10 4 ). As mentioned 
above, the exponential covariance model with quite rough (non-differentiable) realizations and 
the Gaussian model with very smooth (infinitely differentiable) realizations can be obtained 
from the Matern model as special cases for v — 0.5 and v —» oo, respectively. To emphasize the 
flexibility of the correlations developed in the relaxation process we also fit the relatively rough 
data for r = 10 (Fig. 3(a)| ) to the exponential and the very smooth data for r = 10 4 (Fig. 3(d)| ) 
to the Gaussian models. One can see an excellent collapse of both curves with the best fits to 
the Matern model. In spite of the very small value of £ inferred in the Matern model for r = 10 4 
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Figure 4. The objective function O given by © over the £ — v parameter space. O represents 
the fit of the data obtained at T — 10 -4 and r = 10000 (see Fig. 3(d)) to the Matern model 


with a n = 0.001 and a = 0.15. 


(see the explanation below), the characteristic length is expected to increase with decreasing 
temperature and increasing simulation time. This expectation can be justified by the fact that 
the ground state of the MPR model (the equilibrium state at T = 0 corresponding to the lowest 
energy), is the ferromagnetic state with all the spins aligned in the same direction. Thus, for 
T 0 one can expect that the correlation (characteristic) length will span the entire lattice and 
eventually go to infinity for L oo. 

Next, we comment on the big discrepancy between the characteristic length parameters £, 
estimated for the Matern and Gaussian models, for T — 10 -4 and r = 10 4 . In fact, in the Matern 
model the parameter £ has been found to be highly negatively correlated with the smoothness 
parameter v P2U2SU2E]. This correlation is also apparent by looking at the objective function of 
the WLS fit, defined by Eq. ©, which is shown in Fig. [4]in the £ —z/ parameter space, for the data 
presented in Fig. |3(d)| It demonstrates that the samples with small £ and large z> values are quite 
“close” to the samples with large £ and small z> values. Put differently, a small difference between 
samples from the same population (particularly in non-ergodic samples) can lead to completely 
different parameter estimates. Thus, one should be careful when performing parameter inference 
from the empirical variogram. Recently, we have presented a model-independent definition of 
the correlation length borrowed from statistical field theory and proposed to use a so-called 
ergodicity index to compare coarse-grained measures corresponding to both trivial (standard) 
and non-trivial, e.g. Matern, covariance models with different parameters [28]. 

4. Summary and outlook 

We have modified the standard planar rotator spin model from statistical physics to display 
a flexible type of short-range spatial correlations relevant in geostatistical modeling. In 
particular, the empirical study of spin configurations produced by Monte Carlo simulations 
in the nonequilibrium regime at various temperatures and stages shows that the smoothness 
and correlation range vary greatly versus the temperature and the simulation time. This 
behavior implies that the model has good potential for the simulation and prediction of 










spatial processes in geophysical and environmental applications. In particular, one can use 
the model to perform conditional simulations of gridded data, such as remote sensing images. 
Conditional simulations honor the sample values and reconstruct the variability of missing data, 
based on simulation parameters inferred from the available samples. Owing to the fact that 
the model does not show undesirable critical slowing down, the relaxation process is rather 
fast; furthermore, the short-range nature of the interactions between spin variables allows 
vectorization of the algorithm. Consequently, the proposed method is significantly more efficient 
than the conventional geostatistical approaches, and thus applicable to huge datasets, such as 
satellite and radar images. The implementation details are now under investigation. 
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